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o 

^^ ■ We investigate the rheology of sheared granular materials near the jamming transition 

{^JQ' point. We numerically determine the values of the critical fraction and the exponents for the 

^ , jamming transition using a finite size scaling and the nonlinear minimization method known 

.^1^ . as the Levenberg-Marquardt algorithm. The exponents are close to our previous theoretical 

^ ' prediction, but there is a small discrepancy, if the critical point is independently determined. 

o\ , 

'~i. §1. Introduction 

O 

O ■ Athermal disordered materials such as colloidal suspensionsP foamspl' and gran- 

ular materials^ behave as dense liquids when the density is lower than a critical value, 
while they behave as amorphous solids when the density exceeds the critical value. 
This rigidity transition is known as the jamming transition, which could be a key 
. . concept to characterize disorder materials even for glassy materials.'^ 

C^ \ Near the jamming transition point, such materials show critical behavior, where 

the pressure, the elastic moduli, and the characteristic frequency of the density of 
. . state exhibit power law dependences on the distance from the transition point.'^''^''^ 

^ \ In particular, the critical scaling law characterized by a set of critical exponents, sim- 

O ' ilar to those in thermal critical phenomena, is observed in the rheology of athermal 

disordered materials ,® '® '1^ '13 'H^''^^''^"^ ''113 '13 '13 though the transition becomes 
discontinuous under the existence of friction for granular materials.^^ The precise 
values of the critical exponents, however, are still controversial because the values of 



C^ 



O 



(N 
> 

oo 

o 



o 



Q^ '. them are inconsistent among the researchers ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ 

OO I In this paper, we try to numerically determine the critical exponents near the 

jamming transition for granular materials under the plane shear using a nonlinear 
minimization procedure and a finite size scaling for the critical fraction. The contents 
of this paper are organized as follows. Previous results for the critical rheology of 

^Sl ■ athermal disordered materials are summarized in the next section. In § [3l the details 

of our numerical results are presented, where we explain models and their setup in § 
13. H and the critical fraction and the exponents are respectively determined in § 13.21 

K/{ • and ^ 13.31 In ^ 4, we discuss and conclude our results. 

u 

..." §2. Review of scaling properties near the jamming transition 

Let us consider a sheared athermal system characterized by the packing fraction 
(j) and the shear rate 7. We restrict our interest to systems consisting of repulsive 
particles in which the normal interaction force between contacting particles is pro- 
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portional to 6^2 with 612 = r\2 — oyi^ where ryi and a\2 is the distance between the 
particles' center of mass and the average diameter of the particles, respectively. The 
exponent A characterizes the repulsive interaction, i.e. A = 3/2 is for spheres of 
Hertzian contact law, while the simplified linear model {A = 1) is often used. It 
should be noted that the critical properties are determined by the behavior in the 
limit of 5i2 — ?• 0, if the repulsive force cannot be characterized by a single A. Thus, 
when the interaction potential analytic near (5i2 = 0, such a model always belongs to 
the same universality class of Z^ = 1^^ For granular materials, tangential contact 
force exists, but is occasionally ignored to extract universal properties. We call the 
system without the tangential contact force the frictionless system, while the system 
with the tangential force the frictional system. 

It should be noted that the inertia force is always important for granular as- 
semblies, and thus the contact dynamics satisfies an under damped equation. On the 
other hand, the other systems such as foams and colloidal suspensions are believed 
that inertia force is negligible and the contact dynamics is described by an over- 
damped equation. We also note that granular liquids are characterized by Bagnold's 
law in which the pressure P and the shear stress S satisfy 

50C72, Pa7^ (2-1) 

while the other liquids such as dense colloids and foams satisfy Newtonian law 

5oc7, Pcx7. (2-2) 

In the frictionless athermal systems, we believe that the jamming transition is 
continuous. When the packing fraction (j) is lower than the jamming fraction (pj^ 
which is the onset of the rigidity, the system behaves as a liquid. Thus, its rheology 
is characterized by Eq. ()2Tp or ()2-2p depending on the system. When (p is larger than 
(/)j, S and P satisfy 

5 «(</)- (AJ)^^ p oc ((/. - </)J)^^ (2-3) 

with the critical exponents y^ and y',. At the critical fraction (/)j, S and P exhibit 
power laws as 

S<:x^y\ Pcx^y'-i (2-4) 

with the critical exponents y^ and y'. These rheological properties can be rewritten 
as the scaling relationa^J' '11211 '113 



5(7,(/)) =/2'^5 



<p-<Pj 

^/3 



P{lA) = l'''^v{^^y (2-5) 

with the critical exponent /3 = y^/y(j, = y'-^/y'^,- Indeed, to satisfy Eqs. ()2Tp . (|2-2p . 
(|2-3p and ()2-4p . it is sufficient that the scaling functions S{x) and V{x) respectively 
satisfy 

lim S{x) oc x^"*, lim V{x) ex x^'*, (2'6) 
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and 



lim 5(x) oc \x\y^-^/^, lim V{x) oc \xf^-^^^ (2-7) 

re— ^— oo x— >— oo 

for the underdamped system, while 

hm S{x) oc \x\yt'-'^/f^, hm 7'(x) oc |x|^^~^/^ (2-8) 

for the overdamped system. It should be noted that the exponents y,^ and y', are 
believed to be independent of the existence of inertia force. Indeed, the appearance 
of the yield stress is determined only by the force transfer in the percolation network 
of jammed materials. On the other hand, y^ and y' might depend on the detailed 
properties of dynamics. 

Through many simulations and experiments, we recognize that there exist some 
common propevties^^^^^^^^^^^^^^^^ (i) The critical exponents are 
insensitive to the spatial dimension if the dimension is above two, and (ii) the expo- 
nents strongly depend on A. These properties are counter intuitive, and is opposite 
to the conventional critical phenomena. 

Nevertheless, the values of the critical exponents are inconsistent among various 
estimations or observations. In fact, for overdamped frictionless particles with A = 1, 
Olsson and Teitel reported y^ = 1.2 and y^ = 0.42 in their first paper on the 
jamming transition,"^ but in their later paper, ^'^' they estimated y</. = J/i = 1-08 and 
y^ = yiy = 0.28. The theory for overdamped frictionless particles proposed by Tighe 
et ali^ suggests y^ = A + 1/2 and y^ = 1/2, where they assume that the shear 
stress is given by S* = Gjy for <f) > <j)j with the shear modulus^''^ G oc (0 — (/>j) '^ 
and the yield strain 'jy oc (p — (j)j , which give the prediction of y^. Their prediction is 
consistent with the experiment of colloidal suspensions,'^ but contradicts with the 
numerical estimation of Olsson and Teitel.'^ 'l^ 

For the frictionless granular materials with A = 1, the critical exponents are 
reported as y-y = 5/7 and y!y = 4/7 in Ref. 8). Hatano found that the critical 
scaling relation ()2-5p holds with y^ = 1.2, y^ = 0.63, y', = 1.2, and y'^ = 0.57 

in his first report,'^ but y^ and y-y are respectively estimated as 1.5 and 0.6 in 
his recent paper™ Otsuki and Hayakawa proposed a phenomenological theory to 
predict y^p = y'. = A and y^ = y'^ = 2A/{A + 4), but the values differ from Hatano's 

estimation.'^ Note that some of differences between the two groups are superficial. 
Indeed, if we use the same y^, all exponents in one group agree with those of the 
other group. Therefore, the precise estimation of y^ is crucial. 

For frictional granular systems, which are characterized by a microscopic friction 
coefficient ^u, the scaling property (|2-3p using (j)j is no longer valid because the shear 
stress and the pressure change discontinuously at the jamming point. However, by 
introducing a fictitious transition density 1^5 (/u) depending on the friction coefficient 
/i, similar scaling relations exist as 

5(<A,;u)=^(/x){<A-</)5(;u)^^ P{(t^,li)=B{(t>-(Ps{lJ^)Y'^- (2-9) 
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Otsuki and Hayakawa^^' indicate that y^ and y', satisfy y,j, = y's = ^) where the 
prefactor A{ii) depends on /i and B is a constant™ We should note that 4>s{f^) 
coincides (j)j for the frictionless system. 

The estimated values of the exponents in the previous papers are summarized in 
table [H As shown in the table, the values of the exponents differ among the papers. 
Because we expect that the critical exponents y^ and y' characterizing a power law 
liquid depend on the detail of the dynamics, the differences among y^ and y' are 
quite natural. We, however, anticipate that the exponents y^ and y'j^ to characterize 
the quasi static motion are universal. Thus, the discrepancy among the previous 
papers on y^ and y', might be a serious problem. 



Table I. The critical exponents reported in the previous papers, 
system as O, while the underdamped system as U. 



We abbreviate the overdamped 



Paper 
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y<p 


Vi 


y4> 


y-, 


Olsson and Teitel (2007^ 


O (frictionless, A - 


= 1) 


1.2 


0.42 






Olsson and Teitel (201lji^ 


O (frictionless, A - 


= 1) 


1.08 


0.28 


1.08 


0.28 


Tighe, et al. (2010Ji3 


(frictionless) 




Z\+l/2 


1/2 






Nordstrom, et al. (2010J^ 


O (experiment, A — 


3/2) 


2.1 


0.48 






Hatano, Otsuki and Sasa (2007^ 


U (frictionless, A = 


= 1) 




5/7 




4/7 


Hatano (2008J^ 


U (frictionless, A = 


= 1) 


1.2 


0.63 


1.2 


0.57 


Hatano (2010JSJ 


U (frictionless, A = 


= 1) 


1.5 


0.6 






Otsuki and Hayakawa (2009^^*''^^ 


U (frictionless) 




A 


2A 

A+4 


A 


2A 

A+4, 


Otsuki and Hayakawa (2011^^^ 


U (frictional) 




A 




A 





We should note that the estimation of the exponents depends on the choice of the 
critical fraction (^j.HS'lSli Jn Refs. [TT]) . ll5p . ll6p . they simultaneously determined (/)j 
with the critical exponents. However, the critical fraction may have to be determined 
independently as in Refs. UnD,!!!!). In addition, the most of workPS'l^'HS'lIS'l^ 
except for Olsson and TeiteP^J did not use a systematic method, such as the nonlinear 
minimization technique known as the Levenberg-Marquardt algorithm,'^ to estimate 
the critical exponent. 



§3. Numerical result 

Following Olsson and Teitel,'^ we systematically determine the critical expo- 
nents near the jamming point as well as the critical fraction <j)j. In order to determine 
the critical fraction and the exponents, we use a nonlinear minimization technique: 
the Levenberg-Marquardt algorithm. 

3.1. Setup 

Let us consider a two-dimensional granular assembly in a square box with side 
length L. The system includes N grains, each having an identical mass m. The 
position, velocity, and angular velocity of a grain i are respectively denoted by r^, 
Vi, and oji. Our system consists of grains having the diameters 0.7(To, O.Sao, 0.9(To, 
and (To, where there are N/4: for each species of grains. 
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In) 

The contact force f^j consists of the normal part f- and the tangential part 
fij as /jj = fl^ + fl'. The normal contact force /^" between the grain i and 
the grain j is given by f]^' = hl^'O{h]^')0{aij - rij)nij, where hly and riij are 
respectively given by /i^-"' = k^^'{aij — rij) — rj^^'v^'^ and riij = rij/\rij\ with the 
normal elastic constant k^"^', the normal viscous constant rj^"^', the diameter ai of 
grain i, rij = Vi - Vj, cjij = {ai + (Jj)/2 and f^"-^ = {vi - Vj) ■ riij. Here, 0{x) is 
the Heaviside step function defined by 0{x) = 1 for x > and 0{x) = otherwise. 
Similarly, the tangential contact force /^- between grain i and grain j is given 
by the equation /)• = in.m{\h--\,fi\fil\)sign{h--)tij, where min(a,6) selects the 
smaller one between a and 6, and h^- is given by hi- = k^^'u- — rj^^'v- with the 
tangential unit vector tij = {—yij/\rij\,Xij/\rij\). Here, k^^> and r/^*-' are the elastic 
and viscous constants along the tangential direction. The tangential velocity vi- 
and the tangential displacement ui- are respectively given by vi- = {vi — Vj) ■ tij + 
{aiCJi + ajOjj)/2 and ui- = j\^^^f.y^dt viK where "stick" on the integral indicates that 
the integral is performed when the condition \hl-\ < /i|/^" | or another condition 
ufhf} < is satisfied.E3.ED 

We investigate the shear stress S and the pressure P, which are respectively 
given by 

1 /JL_ r A 

(3-1) 






'An) , At)' 

•'ij,y ' ■'ij,y 


\ i j>i 




N 




l^l^^'^r 


f^' + ff^\ 



p=^{y:.i:.r.r{f^'+f^\ , (3-2) 

i j>i I 

where (•) represents the ensemble average. Here, we ignore the kinetic parts of 
S and P, which are respectively given by 5k = — \}2i Pi,xPi,y) /{mV) and Pk = 

X^j Pi' Pi I /(2?7iy), because they are significantly smaller than the potential parts 
in Eqs. (j3-ip and (j3-2p near the jamming transition point. 

In this paper, the shear is imposed along the y direction and macroscopic dis- 
placement only along the x direction by the following two methods. The first method 
is the SLLOD algorithm under the Lees-Edwards boundary condition,'2211 which we 
call "SL" for later discussion, where the time evolution is determined by 

dri Pi , . (o o\ 

-rr = ^ IVif^x, (3-3) 

at m 

^ = E/^J-^P^.^^- (3-4) 

with the peculiar momentum Pj = m{vi — jye^) and the unit vector parallel to the 
x-direction Bx- 
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The second method is quasi-static shearing method, which we call "QS".l2ii'l2^ 
In this method, the shear strain A'j is applied by an affine transformation of the 
position of the particles. Then, the particles are relaxed under the time evolution 
equations 

until the kinetic energy per particles becomes lower than a threshold value Sth- Then, 
we repeat applying the shear and the relaxation process. Here, we chose A^y = 10~^ 
and -Eth = 10~ A;^"'^(Tq, which are small enough not to influence our results. This 
method is expected to correspond to the low shear limit of the SL method. 

In our simulation m, co and A;^"' are set to be unity, and all quantities are 
converted to dimensionless forms, where the unit of time scale is y m-/A;("). We use 
the viscous constants rj^^' = r/'*' = 1.0 and the tangential spring constant k^^' = k^"'' 
for the frictional case. 

3.2. Determination of the critical fraction for the frictionless case 

In this subsection, we determine the transition density <j)j for frictionless particles 
by introducing the jammed fraction / obtained from the simulation using the QS 
method. Here, / is the fraction of samples where the shear stress S is larger than 
a threshold value Sth = 10~^. Figure [J demonstrates the jammed fraction / as a 
function of (p. f is zero in the low density region, while / has a finite value when </> is 
large enough, which suggests the appearance of the yield stress and the rigidity. It is 
to be noted that / around (p = 0.8425 becomes steeper as the system size increases. 
In order to determine (j)j from the data in Fig. [T| we assume f{(j),L) satisfies a 
scaling relation 

/(,/., L) = F((</.-0j)L") (3-7) 

with an exponent a and a scaling function F(x), which satisfies liuix^oo F{x) = 1 and 

\\rax^-ooF{x) = 0. Figure [2] shows the scaling plot based on Eq. (j3-7p . This figure 

confirms the validity of the scaling relation (j3-7|) . Here, we numerically estimate 

(pj = 0.84250 lb 0.00004, where we assume the functional form of the scaling function 

as 

' X + b^ 



F{x) = U + tanh ( — ^ J \ /2 (3-8) 

with the fitting parameters b = 0.0079, Z\x = 0.042, and a = 1.0. Note that the 
critical fraction (j)j = 0.84250 it 0.00004 is almost identical to the simultaneously 
determined value 0.84260 it 0.0004 with the critical exponents.'^''^ However, as 
will be shown, this slight difference between two critical fractions affects the value 
of the critical exponents. 

3.3. Determination of the critical exponents 

In this subsection, let us determine the critical exponents from the simulation of 
the sheared frictionless system using the SL method for 5.0 x 10^^ < 7 < 5.0 x 10~^ 
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Fig. 1. (Color online) Jammed fraction / as a function of <j} ior N ^ 1000, 2000, 4000 and 8000. 
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Fig. 2. (Color online) Scaling plots of the jammed fraction / characterized by Eq. (|3-7|l . The 
solid line is the scaling function given by Eq. (|3-8p . 



with A^ = 4000. Here, we have determined the critical exponents independently, in 
which the critical fraction has been determined as in the previous subsection (case 
A). Figure [3] shows the scaling plots of S and P based on Eq. ()2-5p . This figure 
confirms the validity of the scaling relation (|2-5p . Here, we numerically determine 
y^ = 1.09 ± 0.04, y'^ = 1.06 ± 0.04, and l3 = 0.43 ± 0.01, where we assume the 
functional forms of the scaling functions as 

s(x) = Soil + A,xy^)e{x) + So/ii + Bs\x\^/^-y^)e{-x), (s-g) 

P(2;) = Po(l + Apxy'^)e{x) + Po/(l + Bp\x\'^^^^y'^)e{-x), (3-10) 

which satisfy Eq. ()2-7p with fitting parameters 5*0 = 0.96, Pq = 8.0, Ag = 21 
Ap = 24, Bg = 11157, and Bp = 16803. The estimated values are close to the 
prediction, y^ = 1.0, y'l = 1.0, and /? = 0.4 by Otsuki and Hayakawa,'!^ but a small 
discrepancy exists. It should be noted that Olsson and Teitel reported the exponents 
Ucj} = y's — 1-0^ ^^ -^^f- nHj) . which is close to our results. 

On the other hand, we evaluate the critical exponents y,^ = 1.0 it 0.1, y'f = 
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1.0 ± 0.1, /3 = 0.40 ± 0.01, y-y = 0.40 and y'^ = 0.40 if we simultaneously determine 
both the critical exponents and the critical fraction (case B). We should stress that 
the exponents for case B are identical to those obtained from the mean field theory.l^ 
It is remarkable that the difference of the critical fraction which is about 0.01 % 
affects the values of the critical exponents. We believe that the exponents for case 
A are more appropriate for the critical scaling than those for case B, because the 
critical exponents are only defined in the vicinity of the true critical fraction which 
can be determined independently. In table HIl we compare the exponents for case A 
and case B. 
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Fig. 3. (Color online) (a): Scaling plot of the shear stress S{j,(j)) for the frictionless systems 
characterized by Eq. (|2-5|) . The solid line is the scaling function given by Eq. (|3-9|) . (b): 
Scaling plot of the pressure P(7, <jf>) for the frictionless systems characterized by Eq. (|2-5p . The 
solid line is the scaling function given by Eq, 



3T0|) . Both plots are obtained for case A in 
which the critical fraction is independently determined. 



For the frictional systems, we can only use case B to determine the critical 
exponents based on Eq. (j2-9p from the simulation using the SL method with the 
shear rate 7 = 5.0 x 10^^ and A^ = 4000. This is because the jamming transition 
for frictional grains is discontinuous and the critical exponents are only fictitious 
ones. The estimated values are y^ = 0.97 ± 0.01 and y'. = 0.98 ± 0.01 with the 
fitting parameters (psifi = 0.2) = 0.82, ^^(/i = 0.4) = 0.81, (psifi = 0.8) = 0.79, 
(/>5(/x = 2.0) = 0.78, A{^l = 0.2) = 0.10, A{fi = 0.4) = 0.11, A{n = 0.8) = 0.12, 
A(/i = 2.0) = 0.12, and B = 0.44. The estimated exponents are almost identical to 
those in the previous predictiorP^ and those of frictionless grains for case B. Figure 
m shows the scaling plots of S and P for the frictional particles based on Eq. (|2-9p , 
which verifies the validity of the estimation. 
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Fig. 4. (Color online) (a) : Scaling plot of the shear stress 5(7, (^) characterized by Eq. 
: Scaling plot of the pressure P(7, </!>) characterized by Eq. H2-9|) . 



(b) 



Table II. The critical exponents determined by using a nonlinear minimization method. (Case 
A) : The exponents are determined with r^j obtained in § 13.21 (Case B) : The exponents are 
simultaneously determined with 0,7. 





y</. y'4> P 


Vi Vj' 


<pj 


Frictionless (case A) 
Frictionless (case B) 


1.09 ±0.04 1.06 ±0.04 0.43 ± 0.01 
1.0 ±0.1 1.0 ±0.1 0.40 ±0.01 


0.47 0.46 
0.40 0.40 


0.84250 ± 0.00004 
0.84260 ± 0.0004 


Frictional (case B) 


0.97 ±0.01 0.98 ±0.01 







§4. Discussion and conclusion 



Let us compare our results with those of the previous papers. Tighe et al. pre- 
dicted y^ = 1.5 for the system with A = 1, which is consistent with the numerical 
results for overdamped^ and underdamped systems.'^''^ However, they did use 
any systematic method, such as the nonlinear minimization technique for the de- 
termination of the critical exponents. Our systematic determination of the critical 
exponents for the frictionless granular system gives e.g. y^f) = 1.09 ± 0.04 for cae A 
while case B where the critical exponents are simultanaously determined with the 
critical fraction gives y^ = 1.0 ± 0.1. We should note that y,^ for case A is almost 
identical to another systematic estimation for an overdamped system.'iSJ It still re- 
mains possibility that the estimation depends on the range of the shear rate and the 
density.'^''^ However, our new result for case A may support the suggestion^ that 
y^f) is close but slightly larger than 1. We also note that the previous exponents in 
terms of the mean field theory are almost identical those for case B. It is likely that 
the deviation from the mean field prediction is significant to represent the existence 
of critical fluctuations. 

We should note that the critical scaling of the jamming transition for frictional 
grains is fictitious, because the actual transition is discontinuous. For frictional 
systems, thus, we can only use case B, in which the exponents are almost identical 
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to those for the frictionless case. 

In conclusion, we nuniericahy determined the critical exponents for the jamming 
tranistion of granular materials near the jamming transition point. The estimated 
values for case A are close to the previous theoretical predictiori^ and those for 
case B but a small deviation exists for the frictionless system. The value of case A 
is almost identical to those obtained for the rheology of foams near the transition 
point .1^ The fictitious critical exponents for frictional grains are almost identical to 
those for case B and the theoretical prediction of the frictionless grains. 
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